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Abstract 



Two-dimensional data often have autocovariance functions with elHptical equipotential contours, a property 
known as statistical anisotropy. The anisotropy parameters include the tilt of the ellipse (orientation angle) 
with respect to the coordinate system and the ratio R of the principal correlation lengths. Sample estimates of 
XS ■ anisotropy parameters are needed for defining suitable spatial models and for interpolation of incomplete data. 

The sampling joint probability density fe^R{R^O) characterizes the distribution of anisotropy statistics (R^O). By 
means of analytical calculations, we derive an explicit expression for fe^R{R^O), which is valid for Gaussian, 
stationary and differentiable random fields. Based on it, we derive an approximation /^ r{R-, ^) that is independent 



, . of the autocovariance function and provides conservative confidence regions for the sample-based estimates (i?, 9). 

' ' ' We also formulate a statistical test for isotropy based on the approximation /^ l^{R^6). The proposed /^ l^{R^6) 



provides (i) a stand-alone approximate estimate of the (i?, 0) distribution (ii) computationally efficient initial values 
for maximum likelihood estimation, and (iii) a useful prior for Bayesian anisotropy inference. We validate the 
theoretical analysis by means of simulations, and we illustrate the use of confidence regions with a real-data case 
study. 

Index Terms 



^^ ■ Confidence regions. Machine learning, Monte Carlo simulations. Parameter estimation. Spatial statistics, Stochas- 

: ' tic analysis. 
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CN ■ I. Introduction 



s 



PATIAL random fields (SRFs), also known as spatial random functions, are used in several scientific and 
' engineering disciplines that study spatially distributed processes (e.g., image processing, theory of transport 
^ in heterogeneous media, w^ave propagation in random media, environmental modeling). SRFs w^ith Gaussian joint 
probability density function are also used in machine learning, w^here they are know^n as Gaussian processes |JH. 
For convenience, isotropic SRF models are often used, even though many real data sets display anisotropic patterns. 
Physical anisotropy implies different values of a specific variable along different directions and is expressed by 
means of tensor coefficients (e.g., electrical conductivity tensor). Statistical anisotropy characterizes scalar processes 
(e.g., values of gray-scale images, pollutant concentrations), the correlation range of w^hich depends on the spatial 
direction. Herein wq focus on statistical anisotropy, w^hich implies SRFs w^ith autocovariance functions that possess 
elliptical equipotential contours. Estimation of anisotropy parameters is a topic of ongoing research activity in 
various signal and image processing applications [2J-[9|, as w^ell as in data assimilation. 

In tw^o spatial dimensions, the anisotropy parameters involve the orientation of the principal axes of anisotropy 
and the anisotropy ratio, i.e., the ratio of the principal correlation lengths. A problem of practical interest is 
the estimation of anisotropy parameters from a single available sample. The Covariance Hessian Identity (CHI) 
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method |[TOll , [[TTl is a non-parametric, non-iterative method for obtaining semi-analytic estimates|l| {R^ 9) of SRF 
anisotropy parameters from two-dimensional data sets. Nevertheless, anisotropy estimates are statistics, i.e. random 
variables, the values of which fluctuate between different samples. The fluctuations in the two parameters, as we 
show below, are inter-dependent and non-Gaussian. If the joint distribution of the fluctuations is known, it is possible 
to evaluate whether deviations of anisotropy statistics between different data sets are statistically significant. 

Improved methods for estimating anisotropy are necessary in various research fields. The characterization and 
measurement of anisotropy in biological materials is important for diagnostic and medical reasons [8|, [12]. 
Significant changes in anisotropy over time may suggest a crucial change in the underlying physical processes. For 
example, an accidental release of radioactivity may significantly alter the anisotropy of radioactivity patterns over the 
monitored area. Reliable and computationally fast detection of systematic changes in spatial distributions is crucial, 
especially for automatic monitoring systems [13|. Another practical question is whether an estimated anisotropy 
value actually means significant departure from isotropy; if not, simpler, isotropic autocovariance functions can 
be used for spatial modeling. Answering questions such as the above requires mathematical expressions for the 
confidence regions of anisotropy statistics. Herein we provide closed-form expressions for the confidence regions 
and the probability distribution of anisotropy statistics, which can be used as a prior in Bayesian inference [14]. 

Below, we derive a non-parametric approximation of the sampling joint probability density function (JPDF) of 
anisotropy statistics for differentiable, Gaussian random fields. We prove this expression using CHI, the Central Limit 
Theorem, Jacobi's multivariate transformation theorem, and perturbation expansion. The term "non-parametric" 
implies that the approximation is independent of the SRF autocovariance function (henceforward, covariance 
function for simplicity). The non-parametric approximation is shown to yield a sampling JPDF which is more 
dispersed in parameter space than the exact JPDF. This implies a wider confidence region for the anisotropy 
statistics. Hence, if a sample is classified as isotropic at confidence level p based on the approximate JPDF, it is 
actually isotropic at p' > p. 

The remainder of this manuscript is structured as follows: In Section [III we present essential definitions and an 
overview of CHI, on which the mathematical development of the joint PDF is based. In Section [nil we derive 
a general expression for fo^ji(R^9), using the Central Limit Theorem and the conservation of probability under 
variable transformation. In addition, we obtain an equation that determines the confidence regions of anisotropy 
statistics. In Section[lVlwe derive the non-parametric approximation of fo,R{R^ 0) and the corresponding confidence 
region expression. In Section |V] we formulate a non-parametric test for isotropy. In Section |Vl] we compare the 
non-parametric estimates of confidence regions with the results of numerical simulations, and we illustrate the 
application of confidence regions using real data. Finally, in Section IVIII we review the main results obtained in 
this work, we present our conclusions, and we outline directions for future research. 

II. Preliminaries 

In the following, boldface symbols are used for vectors, matrices and tensors; the superscript "t" denotes the 
transpose of a vector or matrix. Let P C M^ denote the spatial domain and \V\ the enclosed area. The vector s ^V 
denotes the position of a point in V and ||s|| denotes the Euclidean norm of s. Let X(s, uo) represent a scalar SRF on 
the probability space (fi, J^, V) q The state index uj determines the field state and will be suppressed in the following 
for the sake of brevity. The SRF X(s) represents a scalar variable, e.g., temperature, dose rates, or grey-scale 
intensity levels of a digital image. The events in J^ comprise the measured SRF realization(s) or sample states(s). 
E [•] denotes the expectation over the ensemble of states, and the operator Gov (Zi, Z2) = E [Z1Z2] — E [Zi] E [Z2] 
denotes the covariance of the random variables Zi and Z2. 

We will focus on Gaussian SRFs (GSRFs) that possess normal joint probability density functions with variance 
cr^. Wide-sense stationarity will be assumed, i.e. that the mean mx = E [X(s)] is constant, and the covariance 
function Cxx(r) = E [X(s) X(s + r)] — m^ is independent of s. Let e^ denote the unit vector in the spatial direction 
z, (z = 1, 2). If the partial derivatives d'^ c^^{r) / drf in the orthogonal directions z = 1, 2 exist at r = (0, 0), the SRF 



is differentiable in the mean square sense (m.s.s) for every s G P, i.e., lim/^^o E 



\^^X{S) 



h I 



0, 



Vz = 1,2. For Gaussian SRFs, mean square differentiability practically implies that the derivatives of the sample 

^We use the hat over a mathematical symbol to denote random variables that are sample-based estimates (statistics). 
^Q denotes the sample space (ensemble) that includes all the possible states (realizations) of the SRF, T is the set of all observable events, 
J^ C r^, and P{T) G [0, 1] is the probability associated with each event. 
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States exist almost surely ifTSl . We will focus on SRFs that possess at least continuous first-order derivatives in 
the mean square sense. In addition, we will assume short-range correlations, namely that the covariance falls off 
sufficiently fast as ||r|| ^ oo, so that the correlation area Vc'.— -^ J dr c^^{t) is finite. 

Furthermore, we assume that the sample, X^ comprises the values X^ = X(s/e,cj*) of the field X(s) for a 
specific state cj*, where Sj^, k = 1, . . . , A^ are sampling locations. We will denote sample-based estimates by the 
"hat" symbol, e.g., R represents the estimate of the anisotropic ratio R. 

The Covariance Hessian Matrix H (CHM) of a stationary, at least once differentiable, SRF X(s) is defined as 
follows 

^«<"=-w '••'■='■'■ "* 

Let Xij = diX{s) djX{s), z = 1, 2 be the Gradient Kronecker Product (GKP) tensor; X is a symmetric second- 
rank tensor. The expectation of GKP, henceforward called the mean slope tensor and denoted by Q, is defined as 
follows 

Q^J=E[^^X{s)^JX{s)]. (2) 

The mean slope tensor is also known as the matrix of spectral moments, and plays a key role in determining the 
local maxima and excursion sets of random fields [16]. Swerling has proved the following Covariance Hessian 
Identity (CHI) W\' 

Theorem II.l (Swerling's Covariance Hessian Identity). Let X{s) be a statistically stationary GSRF with a 
covariance function that admits partial derivatives d'^ c^^{r) / drf in the orthogonal directions i = 1, 2 a^ r = (0, 0). 
Then, the mean slope tensor is connected to the covariance Hessian matrix as follows: 

Q=H(r)|,^o. (3) 

Based on Jensen's inequality, (E [9iX(s) 92^(s)])^ < E [{diX{s)}'^] E [{52^(s)}^] it follows that det(Q) > 
0, where det(Q) is the determinant of the matrix Q. 

CHI is valid in any number of spatial dimensions. In two dimensions, the anisotropy parameters (i?, 9) satisfy 
the following theorem ifTTIl : 

Theorem II.2. Let X(s) be a statistically stationary and anisotropic GSRF respecting the conditions of \II.l\ and 
^i^i = 1,2 represent the principal correlation lengths of X{s), i.e., ^^ = ^ — ^^r^\r'=o (where a > is an 
0(1) constant). We define the anisotropy ratio R = ^2/^1^ ^^^ ^^^ orientation (rotation) angle, 9, as the angle 
between the horizontal axis of the coordinate system and the first principal axis of the SRF (arbitrarily defined). 
We also define the slope tensor ratios of the elements Qij, i^j = 1, 2.* 

. Q22 l + i^^tan^^ 



Then, the anisotropy parameters are given by 

1 



R = 



Qii 


R^ + tan^e 


Q12 


tan^(i?2_i) 


Qu 


i?2 + tan2 e 


an"^ 


( 2go \ 




Vl-9d/' 


1 


1 -Qd 



9d-(l + 9d)cos2 6i^ 



(5) 



(6) 

n-i/2 

(7) 



Proof. The proof, which is based on Theorem III. 11 is shown in ifTTTl . The only difference is that therein 
R = ^2(1) = C1/C2 is used, while above we defined R = ^2/^1- The results from [llj apply by means of the 
transformation R -^ 1/R. ■ 

Equations (111)-© are invariant under the pair of transformations tan 9 -^ —(tan 9)~^, that is, 6> ^ 6> =b 7r/2, and 
R -^ 1/R. By restricting the parameter space to i? G [0, oc) and 9 G [— 7r/4, 7r/4), or equivalently R G [1, oc) and 
9 G [— 7r/2, vr/2^, the pair {R^ 9) satisfying ©-([5]) for given [q^^ Qo) is unique, thus ensuring that the transformation 
{Qdi Qo) -^ {R^ ^) is one-to-one. Theorem III.2I permits estimating the anisotropy parameters if the mean slope tensor 
can be estimated from the data IfTOll , ifTTll . 



III. Sampling Joint PDF of Anisotropy Statistics 

The sampling values of Qij will be denoted by the random variable Qij. We propose the spatially averaged slope 
tensor estimate Qij, i^j = 1, 2: 

N N 

k=l k=l 

For each sample of the SRF, in general a different estimate of the tensor Qij is obtained, leading through 
application of Q and © to different estimates of R and 9. Hence, a probability distribution is obtained for R 
and 9. In this section we will calculate the joint PDF of R and 9. The calculation of fo^R{R, 9) is based on the 
application of the classical Central Limit Theorem (CLT). As stated below, this implies A^ » 1, a condition that is 
satisfied in most applications of interest. The derivation of fo^R{R, 9) involves several technical steps that include 
the variable transformations from Q to {qaAo) and finally to (R^9). In order to maintain focus in the main text, 
some of the proofs are relegated to Appendices. In Appendix El Jacobi's theorems related to the transformation 
of probability distributions under a variable transformation are reviewed [TSl. In Appendix El we calculate the 
univariate PDFs of the gradient components and the GKP tensor elements. The covariance matrix of the GKP 
tensor elements is derived in Appendix O The JPDF of the GKP tensor elements is derived in Appendix O based 
on the Central Limit Theorem. Appendix |E] presents the proof of Lemma UlI. 31 which determines the bivariate JPDF 
of the average slope tensor ratios, %, and %. Finally, Appendix |Fl proves the explicit relation for the JPDF of the 
anisotropy statistics, ([TSt , as formulated in Theorem IIII.2[ 

A. Central Limit Theorem and joint PDF of 

The classical CLT for scalar random variables is discussed in |[T9ll - l[2T1l . An extension of the classical CLT applies 
to vector random variables [l22l : 

Theorem III.l (Multivariate CLT). Let us assume N independent and identically distributed vector variables 

Zk, k = 1, . . . , A^ with mean m and covariance matrix Czz- The random vector Z = (Zi + • • • + Zn)/N is 
asymptotically (i.e. for N -^ ooj normally distributed with mean m and covariance matrix Czz/N. 

The same theorem also applies to averages of correlated random variables, i.e., random fields Z(s). The main 
requirements are stationarity and finite range of correlations. Stationarity ensures that the correlations between a 
pair of points only depend on the distance but not the location of the pair. If the size of the correlated areas is 
defined as Vc := max f ^^ J dr cz^Zj (r) ) , the requirement is that Vc be finite and \T)\/Vc » 1. This follows by a 

Straightforward extension of the scalar case [1231 . An intuitive explanation of the constraint on the correlation range 
in the scalar case is as follows: consider that V comprises "blobs" of size oc Vc. The point at the center of the blob 
is correlated with other points inside the blob and uncorrelated with points outside the blob. Hence, one can think 
of the blobs as representing independent random variables. Assuming a uniform distribution of sampling points, 
each blob contains, on average, a number of points N^, ^ N Vc/V. Thus, the effective number of "independent 
units" is TVeff = N/N^, ^ \V\/Vc. If K is finite and TVeff > 1, CLT applies to q given by ^. Note that TVeff > 1 
also implies A^ » 1, since A^ = Aeff x A^^. 

B. Joint PDF of slope tensor elements 

We define the vector Q = (Qii? Q22, Qi2)^ that comprises the independent components of the slope tensor 
statistics. According to (|8]), Qij = -^ J2k=i ^iji^k)- The univariate PDFs of the Xij{sk), are derived in Appendix El 

Lemma III.l (Covariance matrix Cg). For a statistically stationary, and anisotropic GSRF, the covariance matrix 

C(g is defined by Cij^kl = Gov [Qij.Qkij for ij, A:, / = 1, 2, i.e., 



Ck 



Cll,ll ^^11,22 ^11,12 
^^22, 11 ^^22, 22 ^^22, 12 
.^12,11 ^12,22 ^12,12. 



(9) 
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The six independent elements (upper triangular entries) of Cg are given by the following series 

= T^ [Qik Qjl + Qi/ Qjk] + ^ 2^ [Hik{Ynm)Hji{Ynm) + Hii{Ynm)Hjk{Tnm)] , (10) 

where Y^r^ = s^ — s^ (n, tti = 1, . . . , A^) /s" ^/ze /ag vector between two locations s^ a^iJ s^. 

Proof: The proof is given in Appendix O ■ 

The leading term ^ [Qi^ Qj7 + Qu Qjk] in (flOl) leads to the non-parametric approximation of fo^R^R^ 9) as shown 
below. The sums over Ynm 7^ contribute parametric corrections that involve the covariance function. By definition, 
Cg is a symmetric matrix, namely Cij^ki = Cki^ij- 

Lemma III.2 (Joint PDF of Q). Let a statistically stationary and anisotropic GSRF X{s) with covariance Cxx(r) 
that is short-ranged and its spectral density satisfies Cxx(k) ^ o{k~^~^) where k = ||k|| and e > 0. Then, the 

k^-oo 

joint PDF of Qu, Q12, and Q22 tends to the following trivariate Gaussian as N ^ 00^ \V\ -^ 00 

where the ensemble mean niQ is given by 

1 ^ 

™^' = ]7^^[(^ll(^^)'^22(Sfc),Xi2(Sfc))] = (Qll,Q22,Ql2), (12) 

k=l 

and the covariance matrix Cg is defined by ^ and (fTOl) . 

Proof: The proof, based on Theorem llll.li is presented in Appendix iDl The condition Cxx(k) ^ o{k~^~^), 

e > implies that V/c ^ oc, 3e > 0, Coo > 0, such that Cxx(k) < Coo/k^^^. This condition is not particularly 
restrictive, since it is satisfied by most finite-range, differentiable covariance functions, including the Gaussian, 
Matem with u > 1 and Spartan covariance models. ■ 

The normal probability plots of Qu, Qu, Q22 of simulated data, shown in Fig. [6l confirm the asymptotic 
normality of the univariate PDFs in agreement with the CLT. 

C. PDF of slope tensor ratios 

Albeit q involves three independent elements, the anisotropic parameters are determined from two slope tensor 
ratios, c.f. ©-([I]). Next, we derive the JPDF of the slope tensor ratios /q(q;niQ,Cg) where q = {qdAoY from 
the joint PDF of Q. 

Lemma III.3 (PDF of slope tensor ratios). For a statistically stationary and anisotropic GSRF X{s) with a 
covariance Cxx(r) that satisfies the conditions of Lemma \III.2\ the PDF /q(q;niQ,Cg) of the slope tensor ratios 
Qd^ Qo is given by the expression 



K _c 

/q(q;mQ,C^) = —^ e 2 



\/2^el^(B2 + 4A) erfcf — ^ | - 4S\/]4 



(13) 



where erfc(-) is the complementary error function, A and B are functions of the random vector q as well as C^ 
and niQ, while C and K are functions of Cg and niQ only: 



^(q,C^)=q*Czlq, (14a) 

B(q,mQ,C^) = -2m*jCziq, (14b) 

C(mQ,C^) = m*jCzimQ, (14c) 

K{C^) = (27r)-3/2 [det(C^)]-V2. (14d) 

Proof: The proof is given in Appendix El ■ 

The joint PDF of q given by (fT3l) is quite different from the Gaussian form of the joint PDF of Q. 

1) Dimensional analysis: Since Qd^Qo are dimensionless, so is /q; the coefficients A^B^C^K, nevertheless, 
have dimensions due to their dependence on Cg and niQ. In addition, ([T3]) impUcitly depends on N through Cg, 
c.f. ([TOl) . Below, we express ([T3l) in terms of dimensionless coefficients and N. 

First, let Cxx(r) = a^p{v;p) where p = (^i,i?, 6>) and — 1 < p(r]p) < 1 is the dimensionless autocorre- 
lation function. Based on ([TOl) and the definition of the covariance Hessian matrix ([T]), it follows that CZ^ = 

A^^^^cr~^P(q;i?, 0), where P(q;i?, 0) is a dimensionless matrix that depends on the functional form of p{r;p). 
Similarly, A^N^^f cj-^ ^(q; R, 9), B ^ N^ il a'^ 0(q; R, 9), C = N^ C(q; R, 9), and K = N^ ^^ a'^ /C(q; i?, 0), 
where the scalar coefficients A^ >B, C, /C depend implicitly on p(r). Based on the above scaling relations, we propose 
the following transformations that involve the dimensionless functions 5(q; i?, 0), (7(q; i?, 6>), ^(q; i?, 6>), which are 
independent of N, a^ and ^i: 

B = 2NV2AB, (15a) 

C = 27V2^, (15b) 

K = V2A^/'^K. (15c) 

Then, the slope-tensor-ratios PDF ([T3l) is expressed as follows: 



/q(q;m(^,C^) = Ke-^^ V^e^^ (25^ iV^ _^ 1) erfc(SiV) - 2S7V 



(16) 



The JPDF and the resulting confidence regions depend on i?, 9 and the functional form of p(r;p). 

2) Asymptotic PDF Limit : Since A > and B < 0, B < and the argument of erfc(S N) in ([T6l) is negative. 
Therefore, for TV ^ oc, and x = 57V, it follows that erfc(x) - 2 + e"^' (-^ + 0(x"2)^ [24J. Hence, to 
leading-order in N, ([T6l) is approximated as follows: 

/q(q;mQ,C^) = 2v^^ (2B^ N' + l) e^^(«^-'5). (17) 

Numerical comparisons show that the absolute relative error between the exact, ([T6]) . and the approximate, (ITtI) . 
JPDF is less than ^ 10~^ for A^ = 50 and ^ 10~^ even for N = 30. The numerical computations are based on 
the non-parametric approximation of Cg derived in Section JVl 

D. Joint PDF of anisotropy statistics 

Theorem III.2 (Joint PDF of anisotropy statistics). For a statistically stationary and anisotropic GSRF X(s) with 
a covariance Cxx(r) that satisfies the conditions of Lemma \IIL2\ the JPDF of the statistics R and 9 is given by the 
following equation: 

op I ^2 1 I 

fe^R{R, 9; mQ, C^) = — -3 /q(q; mQ, C^), (18) 



(r'^cos'^ 9 + sin^' 



where /q(q; niQ, Cg) is given by ([T6l) . 

Proof: The proof is given in Appendix IH 



The function fo,R{R^9) is clearly non-Gaussian; it depends on niQ and Cg via /q(q; niQ, Cg), while q is 
expressed in terms of (R^ 9) using (|4]) and (|5]). If degrees are used, instead of radians, fo,R{R^ 9) should be 
multiplied by tt/ISO, which is the absolute value of the Jacobian determinant of the transformation from radians 
to degrees. 
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Fig. 1: Schematic illustrating the transformation of confidence regions at level p across different coordinate systems. 



E. Confidence Regions 

Herein, we derive expressions for confidence regions of the anisotropy statistics {R^6). Confidence regions 
are used instead of intervals due to the asymmetry of fo^R{R,9). The confidence region for a probability level 
p G [0, 1] is the "volume" of parameter space which contains the sampled values of the statistics with probability 
p. The confidence region is defined by the following equivalent equations: 

7^dQ/^(gii,g22,Qi2;mQ,C^), fcR3 
J^,,dRd9fe^R{R,9;m^,C^). C" c [0,oc) x [-^4,^4) 



The above equations represent the evolution of the confidence region under the variable transformations 
(^, 9) as shown schematically in Fig. [TJ 



q 



Corollary III.l (Parametric equation of confidence regions). For a statistically stationary and anisotropic GSRF 
X(s) with a covariance Cxx(r) that satisfies the conditions of Lemma \IIL2\ the confidence region corresponding to 
level p G [0, 1] in (i?, 9)-space is given by the parametric equation 



B\q; niQ, Cg) - C(q; niQ, Cg) N^ = ln(l - p) 



(19) 



where B, C are given respectively by (I15al) - (ll5bl) and q translates into (i?, 9) by means of & and ([5]). 

Proof: The JPDF /g is given by the trivariate Gaussian (fTTI) . Hence, the confidence region of Q is an ellipsoid 
whose surface satisfies the equation 

)^Czi(( 



nif 



nif 



^P5 



(20) 



where ip = F~^(x^ = p, 2) is the inverse of the chi-square cumulative distribution function with z/ = 2 degrees 
of freedom II251I . Under the transformation Q ^ q, the ellipsoid is projected onto an ellipse, and following the 
transformation q -^ {R^ 9) into an asymmetric convex curve, as shown schematically in Fig. [TJ Based on (IE-461) , 
the equation of the the corresponding ellipsoid in {u^ g^, go)-space is given by 



A(q, Cg) u^ + S(q, mQ, Cg) u + C(mQ, Cg) - . 



0, 



where the coefficients A^B^C slyq given by ([T4]) . For the quadratic equation to have a unique real solution u = Qu 
for any q, the discriminant should vanish, i.e. 



fi2(q;mQ,C(g)-4A(q;C(g) C{mQ,C^)-£p 
We can verify, using ([T4]) . that (1211) represents an ellipse in the space of q, i.e. 



0. 



(21) 



q^ [(C .1 m^) (C .1 m^)^ - (m^ C .^ m^ - £p) C.' 



0. 



Using the dimensionless scaling functions ([T5l) in (|2TI) we obtain the parametric equation 2((7 — B'^)N^ = ^^, 
where by definition F{£p^ z/ = 2) = p; since F(x, z/ = 2) = 1 — exp(— x/2) it follows that £p = —2 ln(l —p) finally 
leading to ([HI). ■ 



IV. Non-parametric JPDF and Confidence Region 

The expressions for the JPDF, fe^R{R, 0; niQ, C^), and the confidence regions of the anisotropy statistics derived 
above depend on the matrix Cg, given by ([TOl) . This matrix involves a series that does not, in general, admit a closed 
form. An approximate, analytical expression can be derived by keeping only the leading term in the covariance 
matrix series (ITOl) . The truncation is justified if Cxx(r) has short-range correlations, which implies that i^zj (r) decays 
fast for ||r|| ^ max(^i,^2)- Then, the leading-order approximation of C^ is given by 



Ck 



2 

N 






12 



Ql2 

Q22 



Q11Q12 
Q12Q22 



_QiiQi2 Q12Q22 2(^12 + ^11^22) 



&^\ 



(22) 



.(0) 



Since C V does not account for the impact of correlations, we expect that it will lead to a joint PDF with higher 
uncertainty, and hence wider confidence regions, than the true PDF. We confirmed this hypothesis by means of 
numerical simulations (see Section IVI-AI) . In addition, C^% as evidenced in (fTOl) . has an 1/A^ scaling prefactor in 

contrast with 1/N'^ for the truncated terms. To accommodate this sample-size dependence, the scaling relations (flSl) 
are accordingly modified below. 

Theorem IV.l (Non-parametric JPDF). For a statistically stationary and anisotropic GSRF X(s) with covariance 
Cxx(r) that satisfies the conditions of Lemma \III.2\ the non-parametric approximation of the JPDF for the anisotropy 
statistics is expressed as 

f(0)/E> n. TD n An _ \Ar.^fl. „M fW( 



where 



f^>^{R, 9- R, e, N) = \deti3g,R)\ /^""{R, 9- R, 9, N), 



f^'\R,9;R,9,N) = 2V^Ko (2BiN + l) eM«o^-V2). 



(23a) 



(23b) 



The coefficients Bq , Kq are given by the following expressions 



and Aq is 



Bo = (2io 



-1/2 



(i?2 - l)(i?2 _ 1) cos(2(^ - 6)) - (i?2 + l)(i^2 + 1^ 



Kq = (vrio)"^/^^^ \{R^ + 1) - {R^ - 1) cos(2^) 



l3 



(24a) 
(24b) 

(24c) 



Aq = {R^ - l)2(i?2 _ 1)2 cos(4(^ - 6)) - 4(^^ - l)(i?^ - 1) cos(2(^ - 6)) 
+ (i?^ + l)(3i?^ + 2R^ + 3) + 2i?2(i?2 _ 1)2 _ 

Proof: Equations (I24al) - (l24cb are derived by replacing Cj| with C^ , defined by (|22)) . in (fT4l) . Also, the 

asymptotic result (ITtI ) of Lemma IIII.3I was used. The Bq and Kq are obtained from the following dimensionless 
scaling functions (I15all - (ll5cl) 



B^2y/NV2ABQ, 

C^2NCo, 
K = V2A^/'^ko. 

In the non-parametric approximation, the coefficient Cq in the exponent is reduced to 1/2 in (I23bb . ■ 

Figures |2] and |3] demonstrate representative plots of the non- parametric JPDF based on (l23T l. Note the bimodal 
structure of the JPDF for A'' = 100 in Fig. |2] with one mode at i? = 1.2 and the other (smaller) at i? «i 0.8. This 
is due to the considerable spread of 9, which results from the relatively small number of sampling points, and the 
degeneracy of the solution, i.e., the fact that the combination {R, 9) is equivalent to (1/i?, 9 — 7r/2); hence, there is 
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(a) N = 100 (b) N = 500 

Fig. 2: Non-parametric JPDF for R = 1.2, 9 = 20° and N = 100, 500 based on (gS]). 




(a) N = 100 




(b) N = 500 



Fig. 3: Non-parametric JPDF for i? = 3, 6> = 10° and N = 100, 500 based on (gJ 



an extended degenerate peak at (0.83, —70°), part of which is folded back into the primary domain. On the other 
hand, for i? = 3 the smaller dispersion of 9 leads to a single mode even for N = 100. 

Corollary IV.l (Non-parametric confidence region). For a statistically stationary and anisotropic GSRF X(s) with 
a covariance c^^(r) that satisfies the conditions of Lemma \III.2\ the non-parametric equation of the confidence 
region corresponding to level p is given by 



o2 1 In(l-P) 



(25) 



2 N 

Proof: The equation is derived by applying the modified scaling relations of Theorem (IIV.II) to the general 
result (EB of Theorem (IIILTI) . ■ 



V. Statistical test of isotropy 

Theorem V.l (Isotropic ratio). Let X(s) be a statistically isotropic GSRF (R = 1) with correlation length ^ 
sampled at N points. In addition, assume that (i) \V\ ^ ^^ and (ii) N ^ 1. The sample values of the statistic R 
are contained with probability p in the following interval: 



Re 



2v^Q^p;Ar(l - ap-N) /l + 2y^ap.N{'^ - Oip-N) 



I -2a 



p]N 



l-2a, 



'p;N 



(26) 



where Ip — F ^(x^ — p^2) — — 21n(l — p) is the inverse of the chi square cumulative distribution function with 
two degrees of freedom and ap.j^ — ^p/N. 



Proof: \V\ » ^^ and A^ » 1 enforce the asymptotic conditions of Corollary lIV.il For i? = 1 the trigonometric 
terms in Bq, i.e., in (I24al) vanish, showing explicitly that the confidence region is independent of 9. Plugging the 
resulting (I24al) in (l25l) , the following parametric equation is obtained: N{R'^ ~ 1)^ ~ 2ip{R^ + 1) = 0. The constraint 
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Fig. 4: Evolution of lower and upper limits of 95% confidence interval for isotropic case {R = 1) versus number 
of sampling points N based on 




Fig. 5: Non-parametric PDF, fR{R), for isotropic case {R = 1) and N = 100 sampling points based on (1231) . 
Shaded area represents the 95% confidence interval, (0.77, 1.29), of R based on (1261) . 



N > 2ip ensures that the roots are real numbers. This is satisfied for A^ » 1; for example, for p = 0.95 it implies 
N > 12. The two admissible (positive) roots of the quartic equation are given by 



Equation (l26l) is independent of the covariance function and thus provides a non-parametric approximation of 
the confidence region for R. The JPDF ([23]) is independent of 9 and ^ for i? = 1. The dependence on N of the 
95% confidence confidence interval for R is shown in Fig. |4l The PDF, fji{R), of i? for i? = 1 and N = 100 
is shown in Fig. [51 including the 95% confidence region predicted by ([26l) . Note that the PDF vanishes, instead 
of peaking, at ^ = 1. This is not an artifact of the non-parametric approximation, since the complete PDF ([TSl) 
also vanishes at ^ = 1. This is due to the root of the Jacobian (IF-471) at ^ = 1, which reflects the fact that the 
single point (1,0) in (g^^, g^) -space is mapped to the straight line ^ = 1 in the (R^ ^) -space. The vanishing of the 
density at ^ = 1 is also borne out in numerical simulations that do not use the Jacobian (see Figure [8d| below). 
This counter-intuitive behavior of fniR), emphasizes the usefulness of Theorem IV. 1 1 as a test for isotropy. 



VI. Application to Simulated and Real Data 

In the following, values of statistics derived from a sample of A^ » 1 points will be denoted by a star superscript, 
e.g. Q*j , ^* , ^* . Average values of the statistics over M different samples will be denoted by a bar over the respective 
symbol, i.e., Q^j. Estimation of Q*. is based on discretized partial derivative operators, dXi(sk), where z = 1,2. 
Discretization introduces errors that increase with the sparsity of the sampling pattern. A "good" sampling pattern 
is characterized by a typical distance a between nearest neighbors which is approximately uniform (ideally, a 
regular lattice pattern is best) and a <^ min(^i,^2), where ^1,^2 are the principal correlation lengths. Different 
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possibilities for dXi(sk) are investigated in [llj. The parameters niQ and C^ in ([TSl) are unknown since they 
represent ensemble properties. For simulated data niQ and Cg are replaced by the averages of the sample estimates, 

J^Q ^ {QiiiQ22iQi2y ^^^ Cg ^ Cq. In the non-parametric approximation, CL^ is obtained from (1221) by 
replacing Qij with Q^j. 

A. Simulated lattice random fields 

We employ multiple SRF realizations with specified (i?, 9) to validate the expression for the confidence region of 
the anisotropy statistics (l25l) . The simulations are conducted on 100 x 100 square grids with lattice constant a = 1, 
by means of the Fourier Filtering Method [26|, f27|. We use Gaussian, Cxx(r) = a^ exp(— r^/^^), and Matem, 
Cxx(r) = cr^ Yi^Y^Ki; (^/O' covariance functions|j. In the Gaussian case, the range of correlations is controlled by 
^ but in Matem case by both ^ and the smoothness parameter v\ the latter adjusts the differentiability of the random 
field between two extremes obtained iox v — \ (non-differentiable exponential function) and v ^ oo (infinitely 
differentiable Gaussian function). For given ^, the field becomes smoother by increasing v. To compensate for this 
effect, we use rescaled correlation lengths ^ = ^/Ic where Ic is the integral scale factor. For Matem correlations 
in d = 2 it follows that Ic — 2\/7n/, while for Gaussian correlations Ic — ^ [28]. 

The non-parametric confidence region obtained by (l25Tl . is compared with the CHI anisotropy estimates obtained 
from 1000 SRF samples. For each sample, we estimate the expectation (|2]) by means of the spatial average (|8]). 
Then, we obtain estimates (^*,^*) of the anisotropy parameters by applying ©-(ItJ. We assume that {R^9) are 
unknown a priori, and we estimate them based on i?, 9. The latter are obtained by calculating the average slope 
tensor, Q over the available realizations and then applying Theorem lll.2[ The use of Q helps to reduce biases due to 
finite grid size and the approximation of derivatives by means of finite differences [11 J. In addition, it compensates 
for potential deviations of the simulated SRF from the target anisotropy values due to the Fourier filtering method. 

In Figs.[6]-[7]we compare the cloud of (^*, 6^*) estimates obtained from each realization with the non-parametric 
confidence region obtained by (l25l) . The latter is denoted by the solid lines that contain the cloud. The estimated 
anisotropy vector, based on Q, is denoted by a small circle inside the cloud. In addition, we include the parametric 
confidence region given by ([T9]) : Cg is estimated by numerical summation of the series (ITOl) using a square window 
function of side ^ 3^ around a grid point s^ at the center of the grid. This restricts the sum over s^o correlated 
neighbors. The summation over all s^ is approximated by multiplying the result with the grid size|j. So long as 
(i) the window area exceeds [2max(^i,^2)]^ and (ii) a < min(^i,^2), the tmncated approximation is stable for 
different box sizes and shapes (e.g., rectangular window). Fig. [6] investigates the isotropic case {R = 1) while the 
anisotropic case i? = 1.5,0 = —30° is considered in Fig. |71 The non-parametric confidence region is more extended 
than the cloud and also encloses the parametric confidence region (contour lines inside the cloud). This result is 
justified, since the non-parametric approximation excludes higher-order correlated terms that reduce the uncertainty 
in Cg. The assumption of slope tensor normality, which was based on the application of CLT in Lemma [ril.2i is 
graphically confirmed by the normal probability plots in Figs. [6el-[6g 



We conducted numerical experiments (not shown here) for several values of the ratio ^/a and various N in 
order to confirm that the non-parametric JPDF is more extended in parameter space more than the actual JPDF. 
When ^/a -^ 0, i.e., as the simulated SRF tends to random uncorrelated noise, the scatter cloud of the anisotropy 
estimates expands and tends to fill the non-parametric confidence region. On the other hand, as ^/a increases, i.e., 
for dense sampling of the SRF, the scatter cloud shrinks inside the corrected confidence region. These observations 
confirm that the non-parametric approximation encloses the actual confidence region. This is expected from an 
information-theoretic viewpoint, since the additional information incorporated in the covariance function should 
lead to less uncertainty (i.e., a tighter confidence region) than the non-parametric approximation which discards 
covariance terms for non-zero distances. 

B. Simulated scattered data 

We simulate scattered data using the following method: First, a realization of an GSRF is generated on a regular 
grid. Then, we randomly choose a fraction of the grid points to mimic scattered data of a real process. Next, 

^ These expressions correspond to the isotropic case. 

"^This approximation introduces errors in the summands from points near the grid boundaries; however, the fraction of boundary points 
for an L X L grid varies as 0{1/L) and is thus negligible for large grids. 
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(a) Scatter plot and CR (Gaussian covariance). 
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(c) Scatter plot and CR (Matem covariance). 
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(b) Random field (Gaussian covariance). 
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(d) Random field (Matern covariance). 
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Fig. 6: (a) and (c): Scatter plot of anisotropy estimates (crosses) and non-parametric confidence region (outer 
contours) for 1000 realizations of isotropic random field with Gaussian (^ = 4) and Matem (^ = 2.51, z^ = 2) 
correlations over an 100 x 100 square grid. A single realization for each field is shown in (b) and (d). Confidence 
regions (CR) are based on anisotropy parameters obtained from the mean slope tensor elements (indicated by small 
circle). Non-parametric confidence regions (outer contours) are obtained by (l25l) . Tighter confidence regions (inner 
contours) are calculated numerically using the appropriate covariance function and incorporate the covariance- 
dependent terms in (fTOl) . Normal probability plots (e)-(g) for the slope tensor estimates (Matem case) justify use 
of CLT in Lemma HIOl 
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(a) Scatter plot and CR (Gaussian covariance) 
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(c) Scatter plot and CR (Matem covariance). 
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(b) Random field (Gaussian covariance). 




(d) Random field (Matern covariance). 



Fig. 7: (a) and (c): Scatter plot of anisotropy estimates and non-parametric confidence region (outer contour) for 
1000 realizations of anisotropic random field with Gaussian (^ = 4) and Matem (| = 2.51, z^ = 2) correlations, 
both having R = 1.5, 9 = —30°. A single realization for each field is shown in (b) and (d). Tighter confidence 
region (inner contour) is numerically calculated and incorporates the covariance-dependent terms in (ITOll . 



we generate several random subsamples from the scattered data set and perform anisotropy estimation for each 
subsample. The subsamples should respect the conditions specified in the opening paragraph of Section |Vll The 
procedure is depicted in Figure |8] Figure ISal demonstrates an GSRF realization on a 1000 x 1000 grid, generated 
by means of the Fourier Filtering Method. The zero-mean, unit- variance isotropic GSRF has a Matem covariance 
with z/ = 2 and ^ = 15. A randomly extracted set of scattered 2000 points is shown in Figure [80 Figure [8cl 
demonstrates a specific subsample of the scattered data set, containing 1000 points. The depicted field corresponds 
to the interpolation of the 1000 points used to estimate (^*,6>*). For the interpolation, we employed the natural 
neighbor interpolation method Il29l of MATLAB®. Details of the impact of different interpolators on anisotropy 
estimation can be found in flTl. Figure [8dl demonstrates the non-parametric confidence region (contour lines) along 
with the scatter cloud of anisotropy estimates, generated by randomly selecting 1000 subsamples of 1000 points 
from the scattered data set. 



C. Environmental emergency scenario 

An application on an environmental emergency scenario data follows. The data represent daily means of radioac- 
tivity gamma dose rates over part of the Federal Republic of Germany, and they were provided by the German 
automatic radioactivity monitoring network [30|. The rates are measured in nanoSievert per hour (nSv/h). The 
normal data set corresponds to typical background radioactivity measurements (^ 100 nSv/h). The emergency data 
includes a simulated local release of radioactivity from the south-west comer of the monitored area that results in 
five stations reporting dose rates around 10 times above the background (exceeding 1000 nSv/h). Table U summarizes 
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(a) Realization of Matern SRF on 1000 x 1000 
square grid. 
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(c) Subsample realization: 1000 randomly selected 
locations (from those in Fig. [8b]) and interpolated 
field (used for anisotropy estimation). 
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(b) Scattered data locations (2000 randomly selected 
points) and interpolated field (for visualization pur- 
poses). 




(d) Anisotropy estimates via subsampling and non- 
parametric CR for N = 1000. 



Fig. 8: Non-parametric confidence region estimation for scattered data, (a) Realization of zero-mean, unit- variance 
isotropic Matern SRF with ^ = 15 and z/ = 2 on a 1000 x 1000 square grid, (b) Locations of random sample that 
includes 2000 points and field realization obtained by means of natural neighbor interpolation, (c) Specific subsample 
realization and interpolated field used for anisotropy estimation, (d) Anisotropy estimates (crosses) generated from 
1000 random subsamples of 1000 points; continuous curve corresponds to 95% non-parametric confidence region 
(CR) calculated with N = 1000 and anisotropy parameters estimated from the mean slope tensor elements as in 
subsection IVI-Al Confidence interval for isotropy is (R^^R^) = (0.925, 1.081). 

TABLE I: Summary statistics of radioactivity dose rate exhaustive data sets (units are in nanoSievert per hour) and 
CHI anisotropy estimates. Abbreviations: min: minimum sample value; max: maximum sample value; std: sample 
standard deviation; sample skewness: sample skewness coefficient; kurtosis: sample excess kurtosis coefficient. 



N = 1008 


min 


mean 


median 


max 


std. dev. 


skewness 


kurtosis 


R* 


r{deg) 


Normal 
Emergency 


57.0 
57.0 


97.7 
106.1 


98.6 
98.9 


180.0 
1528.2 


19.6 
92.5 


0.4 
11.3 


0.6 
144.1 


1.18 
0.45 


7.36 
-0.75 



the statistics of both data sets. The two rightmost columns show the anisotropy parameters estimated with the CHI 
method. The normal set data follows the Gaussian distribution (graph not shown here), and thus has skewness and 
excess kurtosis coefficients close to zero. 

We calculated the JPDF and the 95% confidence regions of the anisotropy statistics based on the anisotropy 
parameters estimated by CHI, i.e., performing natural neighbor interpolation and anisotropy estimation using CHI in 
both exhaustive data sets and then plotting the non-parametric confidence regions and JPDFs. The results are shown 
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Fig. 9: Joint probability density functions (surfaces) and confidence 95% regions (solid contours) for the radioactivity 
dose rate data sets: normal data set (bottom right) and emergency data set (far left). 

in Fig. [9l There is no overlap of the two density functions, and the contours corresponding to the 95% confidence 
regions are non-intersecting, which suggests a statistically significant difference of the anisotropy parameters between 
the background and the emergency data. Since the computation of the anisotropy estimates is very fast, our method 
provides an easy indicator of significant physical change in a system. 



VII. Discussion and Conclusions 

This work focuses on the estimation of anisotropy by means of the Covariance Hessian Identity (CHI) in two- 
dimensional digital data that are either scattered or supported on a grid. We derive explicit expressions for the 
joint PDF of anisotropy statistics, equation (fTSl) , and for the corresponding anisotropy confidence regions at any 
confidence level, equation ([T9]) . The main assumptions used are (i) that the data are drawn from a jointly Gaussian, 
stationary and differentiable random field (ii) that the covariance function is short-ranged. We also derive a non- 
parametric approximation for the joint PDF of the anisotropy statistics, which can be used if the covariance function 
is unknown a priori, or to avoid numerical calculations required for estimating the covariance. The non-parametric 
approximation is given by (l23l) . The corresponding equation for the non-parametric approximation of the confidence 
region is given by (l25l) . Practical application of the results of this research requires the estimation of anisotropy 
statistics using CHI (or other methods). Accurate estimation based on CHI requires a large sample size, N :^ 1 
and a sample domain that is large with respect to the correlation area. The latter may be difficult to satisfy for 
data with large anisotropy (i? » 1 or i? <C 1). In such cases, the CHI estimate tends to underestimate the actual 
anisotropy. We illustrate the application of the joint PDF and the confidence regions with simulated and real data. 
The results of this research can be used to identify significant deviations in anisotropy between data sets, e.g., due to 
structural differences or major changes in the underlying physical process. The computational cost is minimal, since 
the corresponding expressions are analytical. The major computational cost is the estimation of partial derivatives 
of the random field; the user can freely choose a more accurate field derivatives estimator, at the cost of additional 
computational time. Straightforward extension of this work is possible for jointly lognormal data along the lines 
of LlLl. 



Appendix A 
Jacobi's Theorems 

We use Jabobi's theorems [TSl to determine the PDF under the change-of-variables transformations 



q 



Theorem A.l (Jabobi's univariate theorem). Let Z be a continuous random variable and Y — g{Z), where g{Z) 
is a continuous function in Z. If y = g{z) admits at most a countable number of roots Zj = gj (y), j = 1, . . . ,r, 
then 

t V" w -1^ ^^ ^^t^y) 



dy 
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In the case of a multivariate variable transformation, Jacobi's theorem becomes: 

Theorem A.2 (Jabobi's multivariate theorem). Let Z and Y be two n-dimensional continuous random vectors 
with components (Zi, . . . , Z^) and (Yi, . . . , Yn), respectively. The transformation Y = g(Z) represents the set of 
equations yi = QiC^i), I = 1, . . . , n. Assume that the functions gi are continuous and possess continuous partial 
derivatives with respect to each of their arguments. 

(i) If the gi define one-to-one mappings, unique inverse functions g^ such that Z = g~^(Y) exist. Then, the 
transformation of the JPDF fz to /y is accomplished by means of 

fY^fz{g-\Y)) |det(J)|, 

where 3 is the Jacobian of the transformation: 

(ii) IfY — g(Z) admits at most a countable number of roots Tij — gj (y), j = 1, . . . , r, then 

r 

/Y = 5^/z(g7'(y)) |det(J,)|, 
where 3j is the the Jacobian corresponding to the j-th root, defined by 

^ d{yi,...,yn) 

Theorem lA.21 also applies if dim(Y) = m < n. Then, the m-dimensional vector Y is augmented by an (n — m)- 
dimensional vector Z' = h(Z) where h(-) is a simple function with continuous partial derivatives. The n — m 
dummy variables in Z' are eliminated from the JPDF of Y by integration. 

Appendix B 
pdf of gradient product tensor 

Below, we obtain the PDF of the gradient components Xi = diX{s) for z = 1, 2. Let us define by fxXdiX = z) 
the PDF of the gradient component diX, and by fxij{Xij = y) the PDF of the gradient product diX{s)djX{s). 
The following theorem holds [15|: 

Theorem B.l (PDF of field gradient). For a Gaussian, differentiable and stationary SRF X{s), the gradient 
component diX{s) is a zero-mean Gaussian SRF with covariance function given by the following expression: 

E [9,X(s) d,X{s + r)] = -^^^- (B-27) 

In light of ([3]) and (IB-271) the variance of diX{s) is given by Var (5^X(s)) = Qu. Hence, the univariate PDF 
fx, is given by 

fx^ {d,X = z) = —1= e-^V2g. , (B-28) 

A. PDF of diagonal elements of gradient product tensor 

Theorem B.2 (PDF of diagonal elements of gradient product tensor). Let X{s) be a Gaussian, statistically stationary 
SRF that admits (in the mean square sense) all partial derivatives of second-order (at least). Then, the PDF of Xu 
is given by the chi square (xt) distribution with one degree of freedom (z/ = 1).* 

y 

fxAXii^y)^ %^ - (B-29) 



PETRAKIS & HRISTOPULOS: JOINT PDF OF 2D ANISOTROPY STATISTICS 17 

Proof: Let us define Z = 5^X(s), (z = 1, 2) and y = X^^. The equation y = ^(z) = z^ admits two real roots 
for y > 0, i.e., yi^2 = i\/^, but no real roots for y < 0. Since -|— = =^27^' ^^ ^PPlyii^g Theorem lA.il we obtain: 

fxAXii = y) = ^ [/x.(v^) + fxA-^)\ , y > 0. (B-30) 

Equation (IB -291) follows from the above and (IB-281) . The standard density of the Xi distribution is obtained 
from (IB-291) via the replacement y' = y/Qu- ■ 

The mean and variance of Xa are thus obtained by 

E[X,,] = g,„ (B-31) 

\^i{Xu)^2{Quf. (B-32) 

B. PDF of non-diagonal elements of gradient product tensor 

Theorem B.3 (PDF of non-diagonal elements of the gradient product tensor). Let X{s) be an SRF that satisfies 
the conditions of Theorem \B.2\ Then, the PDF of X12 = 5iX(s) 52-^(s) is given by: 

where Q is the slope tensor, det(Q) is its determinant, and Kq is the modified Be ssel function of the second kind 
and order zero. Since Qu = Q21, it holds that /x2i(-^2i — v) — /xi2(-^i2 = y)- 

Proof We define the random vectors Z^ = (aiX(s), a2X(s)) = (Xi,X2), Y^ = (XiX2,X2). The Jacobian 
determinant for the transformation Z ^ Y and its absolute value are given by 



det(J) 



2/2 ^ -yi ^2 ^ 



1 

Thus, according to Theorem IA.2I fx^^ i^ given by 



— ^|det(J)| = ^. 

2/2 1^21 



roc J 

/Xi2(^i2 = 2/1) = / dy2 fzz {yi/y2,y2) 1 — r, (B-34) 

J-00 I2/2I 

where fzz is the bivariate Gaussian PDF of the two-component fluctuation (^1,^2), given by 

Then, the integration in (IB-341) can be evaluated using the integral (3.471.9) LIL P- 368] and the transformation 
yl -^ X, which lead to (IB-331) . ■ 

The mean of X12 is given by (O. The variance is obtained using the normality of the 5iX(s) and 52-^(s) 
distributions by applying the Isserlis-Wick moment factorization theorem ||32]| , Il33l : 

E[Xi2]=gi2, (B-36) 

Var(Xi2) = giig22 + Q?2. (B-37) 

For Ixl -^ oc, Ko(\x\) ^ A/ooe"!^! Ii24il . Hence, based on (IB-33D , fx,. decays for large |Xi2| as 

-IX12I 



^""^^ ^ V2^(gilQ22)V4 y^\X^\ '""^ 



\/Qii Q22 + sign(Xi2) Q12 



(B-38) 



where sign(x) = l,x > A sign(x) = — l,x < 0. 

Asymptotic convergence of fxi2 to zero for |Xi2| -^ 00 requires that the denominator of the exponent be 
positive. Qii and Q22 are positive by definition. Regardless of the sign of Q12, the positive definiteness of Q 
(see Theorem III. IK implies that \/Qii Q22 > ^Qi2- Hence, the denominator is indeed positive. The asymptotic 
dependence of fx^2 matches that of the Xi distribution. 
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Appendix C 
Proof of Lemma HIlT] 

Proof: Using the definition ([8]) we obtain: 

(l ^ 1 ^ \ I 

Cij^kl = Gov I — ^ Xij{Sn), — ^ Xkl{Sm) 1 = ^^ XI ^^^ (X^j (s^), X/e/(s^)) . (C-39) 

\ n=l m=l / n,m 

Due to the translation invariance of X(s), the double series in (IC-39b is reduced to a single series over all (A^^) 
lag vectors r^^ = s^ - s^ (n, m = 1, . . . , TV), i.e., 

Qj,/c/ = j;pYl ^^^ (Aij(so), X/e/(so + r^m)) 

= ^Cov(Xi,-(0),XfcKO)) + ^ 5^ Cov(Xi,(0),Xw(rnm)). (C-40) 

Covariance of the gradient product tensor: Let r denote any lag vector (including r = 0) between two sampling 
points. Based on the definition of the covariance function it follows that 

Gov {Xij{0),Xki{v)) = E [Xij{0)Xki{v)\ - E [Xi,-(0)] E [Xki{v)\ . (C-41) 

Note that E [Xij{0)Xki{r)] = E [^^^(O) djX{0) dkX{r) ^(^(r)]. According to TheoremE] the gradient fields are 
Gaussian SRFs. Hence, E [Xij{0) Xki{r)] can be calculated using the moment factorization property of multivariate 
normal distributions [|321 . [|33l : 

E [XijiO) Xki{r)] = E [diX{0) djX{0)] E [dkX{r) diX{r)] + E [diX{0) dkX{r)] E [djX{0) diX{r)] 

+ E [diXiO) diXir)] E [djXiO) afcX(r)] = i/i,(0) Hki{0) + Hik{r) Hji{v) + Ha{v) Hjk{r). (C-42) 

The last equality follows from the definition of the covariance function. Theorem IB.li and the definition ([T]). The 
second term on the right-hand side of (IC-411) is simply 

E[X,j{0)]E[Xu{r)]=E[d,X{0)djX{0)]E[dkX{T)diX{v)] = H,j{0)Hu^^^ (C-43) 

Thus, in light of (IC-421) and (IC-43L equation (IC-411) becomes 

Cov(X,,(0),X,Kr)) = H,kir)Hji(T) + Hu(r)Hjk(r). (C-44) 

Using (IC-441) in (IC-40L and Theorem III. II to express the zero-lag component of the covariance Hessian matrix, 
equation ([TOl) is obtained. 



Appendix D 
Proof of Theorem IIII. II 

Proof: The Xij{sk) are stationary GSRFs by virtue of the stationarity of the GSRF X{s). Hence, (/)ijki{'^) •= 
Coy{X,j{s),Xki{s + T)) = Cov(X,,(0),Xfc;(r)). Using ([C31, 0.,itKr) = H,k{r)Hji{T) + Ha{r)Hjk{T). The 
range of the GSRF Xij{sk) is determined by the integral 



^^ ^ S (^^ / "^'^^^'^^'^^ I '^^^■'=^(°^ ^ '^) 



Based on (IC-441) . (l)ijki{0) = Qij Q/c/ + Qi/ QjTc and thus (l)ijki{0) has a finite value for finite correlation lengths. 

To calculate /^ dT(l)ijki{r), we assume that |P| ^ oc and express the integral in terms of the Fourier transform 
of Cxx(r). Any permissible covariance function Cxx(r) admits the following pair of transformations, where Cxx(k) 
is the covariance spectral density: 



Cxx(r)=^ydke^'^"ax(k), 



PETRAKIS & HRISTOPULOS: JOINT PDF OF 2D ANISOTROPY STATISTICS 19 

Based on the above, we obtain Hij{v) = j^^ J dkkikj e^^'^C^^(k), and thus 

/ dr(f)ijki{r) = —--2 / dkkikj kkki[C^y,{k)]'^ . 

In the above, j = a/^, k • r = kiri + /c2r2 is the inner vector product, and Jdk = J^^dki J^^dk2 or 
Jdk = J^ kdk Jq^ (i(/) in polar coordinates. Existence of this integral requires integrability of Cxx(k) at /c = 
and at A: ^ oc. Since Cxx(r) is short-ranged, the integral J drc^^{r) = Cxx(O) is finite, and thus the limit fc = 
is well-behaved. At the limit fc ^ oo, the integral converges (using polar coordinates) if [Cxx(k)]^ falls off faster 
than k~^~'^^, where e > 0. This ensures that the covariance (l)ijki{^) is short-ranged. 

Thus, Theorem IlII. II applies to the vector random variable Z^ = ((-^ii(s/c), -^22(8^^)5-^12(8/0))^ leading to (fTTI) . 



Appendix E 
Proof of Lemma [TIO 

Proof: We use Theorem lAl2l with Z ^ Q, Y ^ {Qd^QoY- Since dim(Y) = 2 < dim(Z) = 3, we append to 
Y the dummy variable u = Qu > 0. Using definitions & and ©, the absolute value of the Jacobian determinant 
for the transformation (Qii, Q22, Q12) -^ {qd, Qo, u) is 



Jq = — ^P—- ^ — ^ det(Jq) = u . 

d{qd,qo^u) 

The dummy variable u is integrated according to Theorem IA.2I leading to 

/q(q; mQ, C(g) = / /(g(i/, ^o^, qau] mQ, C(g) i/^ dix. (E-45) 

In terms of q^ and qo, the exponent of the trivariate PDF (fTTI) is transformed as follows 

(Q - niQ)^ C.^ (Q - niQ) = A(q, C^) u^ + S(q, niQ, C^) ix + C(mQ, C^). (E-46) 

By virtue of the above, the integral (IE-451) is expressed as follows (suppressing the dependence of A, S, C, K for 
brevity): 

/q(q; niQ, Cg) = K / u^ ^-\{Au-^Bu^C) ^^ 
Jo 

According to (I14al) . A > because Cg is a covariance matrix, and hence Cg as well as C^^ are positive definiteEl. 
Thus, the Gaussian integral above exists and its value is given by ([T3l) . 



Appendix F 
Proof of Theorem [IIL2] 

Proof: Equation (fTSl) follows from the transformation {q^, qo) -^ {R, 0) with Jacobian matrix Jq^r. According 
to|A21the transformed PDF is given by fo,R(R, 0; mQ, C^) = /q(q; mQ, C^) |det(J^,i^)| where det(J^,i^) is given 
by 

2R (r^ - l) 



det(J^,i^) 



dO OR 
dqo dqo 
dO OR 



3. (F-47) 

^2cos2^ + sin2 0^^ 



Under the restriction of the parameter space to i? G [0, 00) and 9 G [— 7r/4, 7r/4), or equivalently i? G [1, 00) and 
9 G [— 7r/2,7r/2), the transformation {qaAo) -^ (R^^) is one-to-one with a Jacobian determinant given by (IF-47t . 
Finally, based on Theorems IA.2I and lllT31 fe,R{R,9) is given by ([18]). ■ 

^A square p x p matrix M is positive-definite, denoted by M > 0, if x* Mx > V p x 1 vectors x 7^ 0; then, M~^ > 0. 
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